/*==============================================================================
案例 7：信用风险评估（分类算法）
================================================================================

业务场景：
银行需要评估贷款申请人的信用风险，预测是否会违约。
通过分类模型识别高风险客户，降低坏账率，优化信贷决策。

学习目标：
1. 掌握二分类问题的建模流程
2. 学习 Logit、Probit 等分类算法
3. 理解分类模型的评估指标（准确率、精确率、召回率、AUC）
4. 掌握业务阈值的选择

数据来源：nlsw88.dta（模拟为信贷数据）

作者：Stata ML Course
日期：2025-11-03
==============================================================================*/

clear all
set more off
capture log close

* 设置工作目录
cd "`c(pwd)'"

* 创建输出目录
capture mkdir "output/cases"
capture mkdir "output/cases/figures"
capture mkdir "data/cases"

* 开始日志记录
log using "output/cases/case07_credit_risk.log", replace text

display "=========================================="
display "案例 7：信用风险评估"
display "=========================================="
display ""

/*------------------------------------------------------------------------------
第一部分：数据准备（模拟信贷数据）
------------------------------------------------------------------------------*/

display "第一部分：数据准备"
display "------------------"

* 加载数据
sysuse nlsw88, clear

* 删除缺失值
drop if missing(wage) | missing(hours) | missing(ttl_exp) | missing(tenure)

* 创建目标变量：违约标记（基于多个风险因素）
set seed 20251103

* 标准化风险因素
foreach var in wage hours ttl_exp tenure age {
    quietly summarize `var'
    gen `var'_std = (`var' - r(mean)) / r(sd)
}

* 创建违约概率（基于风险评分）
gen risk_score = -0.4 * wage_std ///  // 工资低 -> 风险高
                 -0.2 * ttl_exp_std ///  // 经验少 -> 风险高
                 -0.15 * tenure_std ///  // 任期短 -> 风险高
                 + 0.1 * hours_std ///  // 工作时间长 -> 风险高（不稳定）
                 + runiform() * 0.5  // 随机因素

* 创建二分类目标变量（违约=1，正常=0）
gen default_prob = 1 / (1 + exp(-risk_score))
gen default = (default_prob > 0.3)  // 30% 违约率

label variable default "违约标记（1=违约，0=正常）"

* 创建信贷相关特征
gen monthly_income = wage * hours * 4
gen income_stability = ttl_exp / (age + 1)  // 收入稳定性
gen employment_ratio = tenure / (ttl_exp + 1)  // 就业稳定性
gen high_income = (monthly_income > 2000)
gen stable_job = (tenure >= 3)
gen experienced = (ttl_exp >= 8)

* 创建信用评分（模拟）
gen credit_score = 300 + wage * 10 + ttl_exp * 5 + tenure * 8 + ///
    (married * 50) + (union * 30) + runiform() * 100
replace credit_score = min(credit_score, 850)  // 上限 850

label variable monthly_income "月收入（美元）"
label variable credit_score "信用评分（300-850）"
label variable income_stability "收入稳定性"
label variable employment_ratio "就业稳定性"

display "样本量: " _N
display "违约率: " %5.2f (100 * sum(default) / _N) "%"
display ""

/*------------------------------------------------------------------------------
第二部分：探索性数据分析
------------------------------------------------------------------------------*/

display ""
display "第二部分：探索性数据分析"
display "------------------------"

* 1. 违约分布
tabulate default

graph bar, over(default) ///
    title("违约分布") ///
    ytitle("样本数量") ///
    blabel(bar, format(%9.0f)) ///
    scheme(s2color)
graph export "output/cases/figures/case07_01_default_distribution.png", replace

* 2. 违约率 vs 收入
graph box monthly_income, over(default) ///
    title("月收入分布（按违约状态）") ///
    ytitle("月收入（美元）") ///
    scheme(s2color)
graph export "output/cases/figures/case07_02_income_by_default.png", replace

* 3. 违约率 vs 信用评分
graph box credit_score, over(default) ///
    title("信用评分分布（按违约状态）") ///
    ytitle("信用评分") ///
    scheme(s2color)
graph export "output/cases/figures/case07_03_credit_score_by_default.png", replace

* 4. 不同群体的违约率
bysort married: egen default_rate_married = mean(default)
bysort union: egen default_rate_union = mean(default)
bysort south: egen default_rate_south = mean(default)

display ""
display "不同群体的违约率："
display "-------------------"
tabulate married, summarize(default)
tabulate union, summarize(default)
tabulate south, summarize(default)

* 5. 相关性分析
correlate default monthly_income credit_score ttl_exp tenure age

/*------------------------------------------------------------------------------
第三部分：数据分割
------------------------------------------------------------------------------*/

display ""
display "第三部分：数据分割"
display "------------------"

set seed 20251103
gen random = runiform()
gen train = (random <= 0.75)

count if train
display "训练集样本量: " r(N)
count if !train
display "测试集样本量: " r(N)

* 检查训练集和测试集的违约率
quietly summarize default if train
display "训练集违约率: " %5.2f (r(mean) * 100) "%"
quietly summarize default if !train
display "测试集违约率: " %5.2f (r(mean) * 100) "%"

/*------------------------------------------------------------------------------
第四部分：模型 1 - Logistic 回归
------------------------------------------------------------------------------*/

display ""
display "第四部分：Logistic 回归"
display "-----------------------"

logit default monthly_income credit_score ttl_exp tenure age ///
    married union south if train

estimates store logit_model

* 显示结果
display ""
display "Logistic 回归结果："
display "-------------------"

* 计算边际效应
margins, dydx(*) atmeans

* 预测概率
predict prob_logit if !train, pr

* 预测分类（阈值 0.5）
gen pred_logit = (prob_logit > 0.5) if !train

/*------------------------------------------------------------------------------
第五部分：模型 2 - Probit 回归
------------------------------------------------------------------------------*/

display ""
display "第五部分：Probit 回归"
display "---------------------"

probit default monthly_income credit_score ttl_exp tenure age ///
    married union south if train

estimates store probit_model

* 预测概率
predict prob_probit if !train, pr

* 预测分类
gen pred_probit = (prob_probit > 0.5) if !train

/*------------------------------------------------------------------------------
第六部分：模型 3 - Lasso Logit
------------------------------------------------------------------------------*/

display ""
display "第六部分：Lasso Logit"
display "---------------------"

lasso logit default monthly_income credit_score ttl_exp tenure age grade ///
    married union south industry income_stability employment_ratio ///
    c.monthly_income#c.credit_score c.ttl_exp#c.tenure ///
    if train, selection(cv) rseed(20251103)

* 显示选择的变量
display ""
display "Lasso 选择的变量："
lassocoef, display(coef, standardized)

* 预测概率
predict prob_lasso if !train, pr

* 预测分类
gen pred_lasso = (prob_lasso > 0.5) if !train

/*------------------------------------------------------------------------------
第七部分：模型评估
------------------------------------------------------------------------------*/

display ""
display "第七部分：模型评估"
display "------------------"

* 1. 混淆矩阵和性能指标
foreach model in logit probit lasso {
    display ""
    display "模型: `model'"
    display "============="
    
    * 混淆矩阵
    tabulate default pred_`model' if !train, row
    
    * 计算性能指标
    quietly count if default == 1 & pred_`model' == 1 & !train
    local tp = r(N)
    
    quietly count if default == 0 & pred_`model' == 0 & !train
    local tn = r(N)
    
    quietly count if default == 0 & pred_`model' == 1 & !train
    local fp = r(N)
    
    quietly count if default == 1 & pred_`model' == 0 & !train
    local fn = r(N)
    
    local accuracy = (`tp' + `tn') / (`tp' + `tn' + `fp' + `fn')
    local precision = `tp' / (`tp' + `fp')
    local recall = `tp' / (`tp' + `fn')
    local f1 = 2 * `precision' * `recall' / (`precision' + `recall')
    
    display ""
    display "性能指标："
    display "  准确率 (Accuracy):  " %6.3f `accuracy'
    display "  精确率 (Precision): " %6.3f `precision'
    display "  召回率 (Recall):    " %6.3f `recall'
    display "  F1 分数:            " %6.3f `f1'
}

* 2. ROC 曲线和 AUC
display ""
display "ROC 曲线和 AUC"
display "--------------"

* Logit ROC
quietly estimates restore logit_model
lroc if !train, title("Logistic 回归 ROC 曲线") scheme(s2color)
graph export "output/cases/figures/case07_04_roc_logit.png", replace
local auc_logit = r(area)

* Probit ROC
quietly estimates restore probit_model
lroc if !train, title("Probit 回归 ROC 曲线") scheme(s2color)
graph export "output/cases/figures/case07_05_roc_probit.png", replace
local auc_probit = r(area)

display ""
display "AUC 对比："
display "  Logit:  " %6.4f `auc_logit'
display "  Probit: " %6.4f `auc_probit'

/*------------------------------------------------------------------------------
第八部分：业务阈值优化
------------------------------------------------------------------------------*/

display ""
display "第八部分：业务阈值优化"
display "----------------------"

* 尝试不同的阈值
display ""
display "不同阈值下的性能（Logit 模型）"
display "================================"
display ""
display "阈值    准确率  精确率  召回率  F1分数"
display "----------------------------------------"

foreach threshold in 0.2 0.3 0.4 0.5 0.6 0.7 {
    quietly gen pred_temp = (prob_logit > `threshold') if !train
    
    quietly count if default == 1 & pred_temp == 1 & !train
    local tp = r(N)
    quietly count if default == 0 & pred_temp == 0 & !train
    local tn = r(N)
    quietly count if default == 0 & pred_temp == 1 & !train
    local fp = r(N)
    quietly count if default == 1 & pred_temp == 0 & !train
    local fn = r(N)
    
    local acc = (`tp' + `tn') / (`tp' + `tn' + `fp' + `fn')
    local prec = `tp' / (`tp' + `fp')
    local rec = `tp' / (`tp' + `fn')
    local f1 = 2 * `prec' * `rec' / (`prec' + `rec')
    
    display %4.2f `threshold' "    " %6.3f `acc' "  " %6.3f `prec' "  " %6.3f `rec' "  " %6.3f `f1'
    
    drop pred_temp
}

/*------------------------------------------------------------------------------
第九部分：业务应用
------------------------------------------------------------------------------*/

display ""
display "第九部分：业务应用"
display "------------------"

* 1. 风险分层
display ""
display "1. 客户风险分层"
display "----------------"

gen risk_level = .
replace risk_level = 1 if prob_logit <= 0.2 & !train  // 低风险
replace risk_level = 2 if prob_logit > 0.2 & prob_logit <= 0.4 & !train  // 中低风险
replace risk_level = 3 if prob_logit > 0.4 & prob_logit <= 0.6 & !train  // 中高风险
replace risk_level = 4 if prob_logit > 0.6 & !train  // 高风险

label define risk_lbl 1 "低风险" 2 "中低风险" 3 "中高风险" 4 "高风险"
label values risk_level risk_lbl

tabulate risk_level if !train

* 各风险层级的实际违约率
display ""
display "各风险层级的实际违约率："
tabulate risk_level default if !train, row

* 2. 信贷决策建议
display ""
display "2. 信贷决策建议"
display "----------------"

gen credit_decision = ""
replace credit_decision = "批准" if risk_level == 1 & !train
replace credit_decision = "批准（降额）" if risk_level == 2 & !train
replace credit_decision = "人工审核" if risk_level == 3 & !train
replace credit_decision = "拒绝" if risk_level == 4 & !train

tabulate credit_decision if !train

* 3. 预期损失估算
display ""
display "3. 预期损失估算"
display "----------------"

* 假设：贷款金额 10 万元，违约损失率 60%
local loan_amount = 100000
local loss_rate = 0.6

gen expected_loss = prob_logit * `loan_amount' * `loss_rate' if !train

summarize expected_loss if !train, detail

display ""
display "预期损失统计（元）："
display "  平均: " %10.0f r(mean)
display "  中位数: " %10.0f r(p50)
display "  最大: " %10.0f r(max)

/*------------------------------------------------------------------------------
第十部分：保存结果
------------------------------------------------------------------------------*/

display ""
display "第十部分：保存结果"
display "------------------"

* 保存预测结果
preserve
keep if !train
keep default prob_logit prob_probit prob_lasso pred_logit risk_level ///
    credit_decision expected_loss monthly_income credit_score
export delimited using "output/cases/case07_predictions.csv", replace
display "预测结果已保存: output/cases/case07_predictions.csv"
restore

/*------------------------------------------------------------------------------
总结
------------------------------------------------------------------------------*/

display ""
display "=========================================="
display "案例 7 完成！"
display "=========================================="
display ""
display "主要成果："
display "1. 建立了 3 个信用风险分类模型"
display "2. Logit 模型 AUC = " %6.4f `auc_logit'
display "3. 完成风险分层和信贷决策建议"
display "4. 提供预期损失估算"
display "5. 生成 5 个可视化图表"
display ""
display "关键发现："
display "- 收入和信用评分是违约的主要预测因素"
display "- 最优阈值可根据业务目标调整"
display "- 风险分层可支持差异化信贷策略"
display ""
display "输出文件："
display "- 图表: output/cases/figures/case07_*.png (5个)"
display "- 预测: output/cases/case07_predictions.csv"
display "- 日志: output/cases/case07_credit_risk.log"
display ""

log close

